
###################################
### Comparative Party Manifesto ###
###################################

rm(list = ls())
gc()

library(tidyverse)
library(stargazer)
library(haven)
library(data.table)
library(lfe)
library(BBmisc)
library(msm)

### Set WD
setwd("~/Dropbox (Princeton)/Research/*Spatial_Inequality/paper_SMD_inequality/JoP_final/replication")  


### Load CMP data
load("data/cmp.RData")

### CPSD
load("data/cpds.RData")

cpds <- cpds[ ,.(ISO3,country,year,rae_leg)]

d_cmp <- merge(cmp, cpds, 
               by.x=c("cty_code","year"),
               by.y=c("ISO3","year"), all.x=T)

### WB
load("data/wb.RData")

d_cmp <- merge(d_cmp[ year >= 2000 ], wb, by=c("cty_code","year"), all.x=T)

### IDEAS Turnout
load("data/turnout.RData")

### Take turnout average in countries with both parliamentary and presidential elections
turnout <- turnout[ , lapply( .SD, function (x) mean(x, na.rm=F) ), by=c("countryname","year"), .SDcols="turnout" ]

d_cmp <- merge(d_cmp, turnout, by=c("countryname","year"), all.x=T)

### Regional income
load("data/inc_comb.RData")

### SD
inc_comb_sd <- inc_comb[ , lapply(.SD, function(x) sd(x, na.rm=T) ) , by=c("countryname","year"), .SDcols=c("income") ]
setnames(inc_comb_sd, "income","sd_income")

inc_comb_sd <- inc_comb_sd[ !is.na(sd_income) & countryname %in% unique(d_cmp$countryname) ]

setkey(inc_comb_sd, countryname, year)
inc_comb_sd[ , year := as.character(year) ]

# Calculate percentage change from one year to another
inc_comb_sd[ , sd_income_lag_1yr := shift(.SD, type='lag', n=1), by = "countryname", .SDcols="sd_income"]
inc_comb_sd[ , sd_income_pct_chg := (sd_income - sd_income_lag_1yr) / sd_income_lag_1yr ] 

### Merge
d_cmp <- merge(d_cmp, inc_comb_sd, by=c("countryname","year"), all.x=T)


### ---------------
### Prepare data 
### ---------------

### Code party type
d_cmp[ parfam %in% c(10,20,30) , party_type_v2 := "left-wing" ] # includes ecological parties
d_cmp[ parfam %in% c(40) , party_type_v2 := "liberal" ]
d_cmp[ parfam %in% c(50,60,70) , party_type_v2 := "right-wing" ] # with nationalist
d_cmp[ , party_type_v2 := factor(party_type_v2, levels=c("right-wing","liberal","left-wing"))]

### Electoral rule
d_cmp[ , plurality := ifelse(countryname %in% c("France","United Kingdom","United States","Japan","Lithuania"), 1, 0) ]
d_cmp[ countryname == "Italy" & year %in% c(1993:2004), plurality := 1 ]

d_cmp[ countryname == "Hungary" & year <= 2011, plurality := 0 ]
d_cmp[ countryname == "Hungary" & year >= 2012, plurality := 1 ]

d_cmp[ , seat_share := absseat / totseats ]
d_cmp[ , pervote := pervote / 100 ]

d_cmp[ , welfare_norm := BBmisc::normalize(welfare, method="range", range=c(0,1)) ]

d_cmp[ , gdp_log := log(gdp) ]


##########################
### Summary statistics ###
##########################

d_cmp_sub <- na.omit(d_cmp[ ,.(countryname, year, partyname, welfare, welfare_norm, party_type_v2, plurality, 
                               sd_income, pervote, seat_share, turnout, rae_leg, gdp_log, unempl_rate) ])

sum_stat <- as.data.table(rbind(cbind("var" = "Welfare position", "mean" = mean(d_cmp_sub$welfare), "sd" = sd(d_cmp_sub$welfare), "min" = min(d_cmp_sub$welfare), "max" = max(d_cmp_sub$welfare)), 
                                cbind("var" = "Party vote share", "mean" = mean(d_cmp_sub$pervote), "sd" = sd(d_cmp_sub$pervote), "min" = min(d_cmp_sub$pervote), "max" = max(d_cmp_sub$pervote)),
                                cbind("var" = "Party seat share", "mean" = mean(d_cmp_sub$seat_share), "sd" = sd(d_cmp_sub$seat_share), "min" = min(d_cmp_sub$seat_share), "max" = max(d_cmp_sub$seat_share)),
                                cbind("var" = "SD regional income", "mean" = mean(d_cmp_sub$sd_income), "sd" = sd(d_cmp_sub$sd_income), "min" = min(d_cmp_sub$sd_income), "max" = max(d_cmp_sub$sd_income)),
                                cbind("var" = "Election turnout", "mean" = mean(d_cmp_sub$turnout), "sd" = sd(d_cmp_sub$turnout), "min" = min(d_cmp_sub$turnout), "max" = max(d_cmp_sub$turnout)),
                                cbind("var" = "Legislative fractionalization index", "mean" = mean(d_cmp_sub$rae_leg), "sd" = sd(d_cmp_sub$rae_leg), "min" = min(d_cmp_sub$rae_leg), "max" = max(d_cmp_sub$rae_leg)),
                                cbind("var" = "GDP (log)", "mean" = mean(d_cmp_sub$gdp_log), "sd" = sd(d_cmp_sub$gdp_log), "min" = min(d_cmp_sub$gdp_log), "max" = max(d_cmp_sub$gdp_log)),
                                cbind("var" = "Unemployment rate", "mean" = mean(d_cmp_sub$unempl_rate), "sd" = sd(d_cmp_sub$unempl_rate), "min" = min(d_cmp_sub$unempl_rate), "max" = max(d_cmp_sub$unempl_rate))))

sum_stat[ , 2:5 := lapply(.SD, function(x) round(as.numeric(as.character(x)),2) ), .SDcols= 2:5 ]


### SI Table E.1
print(xtable::xtable(sum_stat, include.colnames=T, caption="Summary statistics"), include.rownames=F)



### -------------
### Models
### -------------

d_cmp_lw <- unique(na.omit(d_cmp[ party_type_v2 == "left-wing", .(welfare_norm, plurality, sd_income_pct_chg, unempl_rate, pervote, seat_share, 
                                                                  turnout, rae_leg, year, countryname, partyname, edate, gdp_log) ]))

m1.1 <- felm(welfare_norm ~ plurality*sd_income_pct_chg | countryname + year | 0 | countryname, data = d_cmp_lw)
m1.2 <- felm(welfare_norm ~ plurality*sd_income_pct_chg + gdp_log + unempl_rate | countryname + year | 0 | countryname, data = d_cmp_lw)
m1.3 <- felm(welfare_norm ~ plurality*sd_income_pct_chg + gdp_log + unempl_rate + pervote + seat_share + turnout + rae_leg | countryname + year | 0 | countryname, data = d_cmp_lw)
m1.4 <- felm(welfare_norm ~ plurality*sd_income_pct_chg + gdp_log + unempl_rate + pervote + seat_share + turnout + rae_leg | countryname + partyname + year | 0 | countryname, data = d_cmp_lw)

covar.lab1 <- c("Plurality rule","Change in regional inequality",
                "GDP (log)","Unemployment rate","Party vote share","Party seat share",
                "Turnout","Legislative fractionalization index",
                "Plurality rule $\\times$ Change in regional inequality")

### Table 4
stargazer(m1.1, m1.2, m1.3, m1.4,  
          digits.extra=0, digits=2, column.sep.width = "-10pt", 
          align=T, no.space=T, omit.stat=c("ser","F"), type="latex",
          dep.var.labels = c("Welfare position"), 
          keep = c("plurality","sd_income_pct_chg"),
          covariate.labels = c("Plurality rule","Change in regional inequality", 
                               "Plurality rule $\\times$ Change in regional inequality"),
          title = "Effect of Change in Regional Inequality on Left Party Welfare Position by Electoral Regime",
          table.placement = "h!", label = "tab:manifesto_mod1x",
          add.lines = list(c("Mean DV",
                             round(mean(m1.1$fitted.values),2),
                             round(mean(m1.2$fitted.values),2),
                             round(mean(m1.3$fitted.values),2),
                             round(mean(m1.4$fitted.values),2)),
                           c("Country FE","$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                           c("Election year FE","$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                           c("Party FE","-","-","-","$\\checkmark$")))

### SI Table E.2
stargazer(m1.1, m1.2, m1.3, m1.4,  
          digits.extra=0, digits=2, column.sep.width = "-10pt", 
          align=T, no.space=T, omit.stat=c("ser","F"), type="latex",
          dep.var.labels = c("Welfare position"), 
          covariate.labels = covar.lab1,
          title = "Effect of Change in Regional Inequality on Left Party Welfare Position by Electoral Regime",
          table.placement = "h!", label = "tab:manifesto_mod1x",
          add.lines = list(c("Mean DV",
                             round(mean(m1.1$fitted.values),2),
                             round(mean(m1.2$fitted.values),2),
                             round(mean(m1.3$fitted.values),2),
                             round(mean(m1.4$fitted.values),2)),
                           c("Country FE","$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                           c("Election year FE","$\\checkmark$","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                           c("Party FE","-","-","-","$\\checkmark$"),
                           c("Economic covariates","-","$\\checkmark$","$\\checkmark$","$\\checkmark$"),
                           c("Political covariates","-","-","$\\checkmark$","$\\checkmark$")))





### Plot
estmean.m1.3 <- m1.3$coefficients
var.m1.3 <- m1.3$clustervcv

rownames(estmean.m1.3)

marg.eff1 <- estmean.m1.3[2]
marg.eff2 <- estmean.m1.3[2] + estmean.m1.3[9]

se1 <- msm::deltamethod (~  (x2), estmean.m1.3, var.m1.3, ses=T)
se2 <- msm::deltamethod (~  (x2) + (x9), estmean.m1.3, var.m1.3, ses=T)

out.m1.3 <- as.data.table(rbind(cbind(coef=marg.eff1, se=se1, type="Proportional\nrepresentation"),
                                cbind(coef=marg.eff2, se=se2, type="Plurality\nrule")))

out.m1.3[ , 1:2 := lapply(.SD, function(x) as.numeric(x)), .SDcols= 1:2 ]

out.m1.3[ , lower90 := coef-1.645*se]
out.m1.3[ , lower95 := coef-1.96*se]
out.m1.3[ , upper90 := coef+1.645*se]
out.m1.3[ , upper95 := coef+1.96*se]


### Figure 8
pdf("out/fig8.pdf", width=4.5, height=2.5)  
(plot <- (ggplot(data=out.m1.3, aes(y=coef, x=type, colour=type, fill=type, shape=type))
          + geom_hline(yintercept=0, colour = "gray50", linetype="dashed", linewidth=.5)
          + geom_errorbar(aes(ymin=lower90, ymax=upper90), width=0, size=1.25, show.legend = F)  
          + geom_errorbar(aes(ymin=lower95, ymax=upper95), width=0, size=0.5, show.legend = F)
          + geom_point(size = 3.5, colour="black", show.legend = F)
          + theme_bw()
          + coord_flip()
          + ylab("Marginal effect of change in\nregional inequality on party welfare position")
          + xlab("Electoral regime")
          + theme(legend.position="right", panel.border = element_rect(color="gray"), 
                  panel.grid = element_blank(), axis.title.y = element_blank())
          + scale_colour_manual(values=c("red","blue","black"))
          + scale_fill_manual(values=c("red","blue","black"))
          + scale_shape_manual(values=c(21,22,23)) 
))
dev.off()





